1 Carga de paquetes

Son muchos los paquetes empleados en estos análisis. Puedes consultar en el ChatGPT qué hace cada uno. Considera un aspecto también importante: algunas funciones escritas por mí se cargan con source_url y source; dentro de algunas de dichas funciones, también se cargan paquetes adicionales.

library(vegan)
library(sf)
library(tidyverse)
library(tmap)
library(kableExtra)
library(broom)
library(cluster)
library(gclus)
library(pvclust)
library(foreach)
library(leaps)
library(caret)
library(RColorBrewer)
library(indicspecies)
library(dendextend)
library(adespatial)
library(SpadeR)
library(iNEXT)
library(GGally)
library(vegetarian)
library(leaflet)
library(leaflet.extras)
library(readxl)
r <- 'R/'
gh_content <- 'https://raw.githubusercontent.com/'
gh_zonal_stats <- paste0(gh_content,
                         'geofis/zonal-statistics/0b2e95aaee87bf326cf132d28f4bd15220bb4ec7/out/')
repo_analisis <- 'biogeografia-master/scripts-de-analisis-BCI/master'
repo_sem202202 <- 'biogeografia-202202/material-de-apoyo/master/practicas/'
repo_sem202302 <- 'biogeografia-202302/practicas/main/'
devtools::source_url(paste0(gh_content, repo_analisis, '/biodata/funciones.R'))
devtools::source_url(paste0(gh_content, repo_sem202202, 'train.R'))
devtools::source_url(paste0(gh_content, repo_sem202202, 'funciones.R'))
fuentes_manuscrito <- 'fuentes/manuscrito/'
source(paste0(gh_content, repo_sem202302, r, 'funciones.R'))
umbral_alfa <- 0.05

2 Leer y preparar datos

2.1 Formularios de Bryan Josué Fúnez Ramos

nombre_odk <- estudiantes_todos_datos %>% 
  filter(grepl(params$estudiante, `Nombres y apellidos compatible Params`)) %>%
  pull(odk_user)
mis_forms_campo <- odk_campo_form_usuario %>% 
  filter(grepl(nombre_odk, usuario_ODK)) %>% 
  inner_join(odk_campo) %>% 
  select(KEY, usuario_ODK, hexagono, Latitude, Longitude,
         Altitude, Accuracy, fechahora, llovio_ultima_hora,
         foto, responsable, otra_persona, observaciones_finales)
mis_forms_id <- odk_id_form_usuario %>% 
  filter(usuario_ODK == nombre_odk) %>% 
  inner_join(odk_id) %>% 
  select(-SubmissionDate, -`meta-instanceID`)

2.2 Mapa de mis puntos colectados

Nota: sólo figuran los puntos con coordenadas. El total de formularios podría ser mayor al número de puntos mostrado en el mapa. Para un mapa comprensivo, ver más adelante el mapa bajo el texto “Mapa de hexágonos visitados”.

tryCatch(
  expr = mis_forms_campo_sf <- mis_forms_campo %>%
  filter(!is.na(Latitude) | !is.na(Longitude)) %>% 
  st_as_sf(coords = c('Longitude', 'Latitude')),
  error = function(e) message('Se produjeron uno o varios errores al generar el objeto espacial. Esto suele deberse a que a alguno o a todos los puntos de campo les faltaron coordenadas. Es probable que el mapa no se genere o que le falten puntos.'),
  warning = function(warn) {
  message('Se produjeron una o varias advertencias al generar el objeto espacial. Esto suele deberse a que a alguno o a todos los puntos de campo les faltaron coordenadas. Es probable que el mapa no se genere o que le falten puntos.', appendLF = TRUE)  
})
tryCatch(
  expr = leaflet(mis_forms_campo_sf) %>%
    addCircleMarkers(
      radius = 8, popup = ~ paste0(hexagono), stroke = T,
      weight = 1, fillColor = 'white', color = 'black', fillOpacity = 1,
      label = ~ hexagono,
      labelOptions = labelOptions(
        noHide = TRUE, direction = 'auto', offset = c(10, 0), textOnly = T,
        style = list('color' = 'black', 'font-weight' = 'bold',
                     'font-size' = '14px'))) %>%
    addTiles(group = 'OSM') %>%
    addProviderTiles("Esri.NatGeoWorldMap", group="ESRI Mapa") %>%
    addProviderTiles("Esri.WorldImagery", group="ESRI Imagen") %>%
    addProviderTiles("CartoDB.Positron", group= "CartoDB") %>%
    addLayersControl(
      baseGroups = c("CartoDB", "ESRI Imagen", "OSM", "ESRI Mapa"),
      position = 'bottomright',
      options = layersControlOptions(collapsed = FALSE)) %>% 
    addFullscreenControl(),
  error = function(e) message('Se produjeron uno o varios errores al generar el mapa. Esto suele deberse a que a alguno o a todos los puntos de campo les faltaron coordenadas. Es probable que el mapa no se haya generado o que le falten puntos.'),
  warning = function(warn) {
  message('Se produjeron una o varias advertencias al generar el mapa. Esto suele deberse a que a alguno o a todos los puntos de campo les faltaron coordenadas. Es probable que el mapa no se haya generado o que le falten puntos.', appendLF = TRUE)  
})

2.3 Fotos

1696777317951.jpg

1696779427975.jpg

1696782836660.jpg

1696789483838.jpg

1696819374386.jpg

1696961966001.jpg

1696962513950.jpg

1696965860056.jpg

1697043574501.jpg

1697048065094.jpg

1697056747040.jpg

1697059774961.jpg

1697834159724.jpg

2.4 Generar matriz de comunidad y ambiental, prepararlas para cargarlas

3 Análisis exploratorio de datos (AED)

3.1 Cargar la matriz de comunidad

spp_odk <- read_excel(
  path = 'odk/biogeografia_hormigas_202302_identificacion.xlsx', sheet = 2) %>% 
  select(list_name, name, label) %>% 
  filter(list_name == 'especieid')
mc <- mis_forms_id %>%
  select(codigo_hexagono_elegido, codigo_hexagono_otro,
         especieid, especieidotra) %>% 
  mutate(hexagono = coalesce(codigo_hexagono_elegido, codigo_hexagono_otro)) %>% 
  select(-matches('codigo_hexagono_.*')) %>% 
  mutate(especieid_todas = case_when(
    !is.na(especieidotra) ~ paste(especieid, especieidotra, sep = " "),
    TRUE ~ as.character(especieid))) %>% 
  separate_rows(especieid) %>%
  distinct() %>% 
  left_join(spp_odk, by = c("especieid" = "name")) %>% 
  select(hexagono, especie = label) %>% 
  filter(! especie %in% c('reina(s)', 'Otra')) %>%
  mutate(presencia = 1) %>% 
  arrange(especie) %>% 
  distinct() %>% 
  pivot_wider(names_from = especie, values_from = presencia, values_fill = 0) %>% 
  arrange(hexagono)
write.csv(mc,
          paste0(fuentes_manuscrito, 'matriz-comunidad-', params$estudiante, '.csv'),
          row.names = F)
mc <- read.csv(
  file = paste0(fuentes_manuscrito, 'matriz-comunidad-', params$estudiante, '.csv'),
  row.names = 'hexagono', check.names = F)
mc %>% estilo_kable(
  titulo = 'Matriz de comunidad',
  nombres_filas = T, alinear = 'r')
TABLA 3.1: Matriz de comunidad
Brachymyrmex heeri Brachymyrmex obscurior Cardiocondyla emeryi Cardiocondyla minutior Cyphomyrmex rimosus Cyphomyrmex sp. Dorymyrmex antillana Monomorium floricola Monomorium pharaonis Odontomachus bauri Paratrechina longicornis Pheidole subarmata Pseudomyrmex subater Solenopsis geminata Tapinoma litorale Tapinoma melanocephalum Tetramorium bicarinatum Tetramorium caldarium Tetramorium lanuginosum
H0268 1 0 0 0 0 0 1 0 0 0 0 1 1 1 0 0 0 0 0
H0517 1 0 1 0 0 1 1 1 0 0 1 0 0 0 0 0 0 0 1
H0571 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0
H0647 1 0 0 0 0 0 1 1 1 0 0 0 0 1 0 0 1 0 0
H0736 1 1 0 0 0 0 1 0 1 0 1 1 0 0 0 0 0 0 1
H0834 1 0 0 0 1 0 0 0 0 1 0 1 0 1 0 0 0 0 1
H0880 1 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 1
H0911 1 0 1 0 1 0 1 0 0 0 1 1 1 1 0 1 0 0 1
H0944 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1
H0956 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0
H1069 1 0 1 1 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0
H1119 1 1 0 0 0 0 1 1 0 0 1 0 0 0 0 0 1 0 1
H1207 1 0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 1 0 0
H1304 0 0 0 0 0 0 1 0 0 0 1 1 0 1 0 0 0 1 1
H1344 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 1 0 0
data.frame(Especies = sort(names(mc))) %>%
  estilo_kable(titulo = 'Lista de especies', cubre_anchura = F, alinear = 'c') %>% 
  column_spec(column = 1, width = "15em")
TABLA 3.2: Lista de especies
Especies
Brachymyrmex heeri
Brachymyrmex obscurior
Cardiocondyla emeryi
Cardiocondyla minutior
Cyphomyrmex rimosus
Cyphomyrmex sp.
Dorymyrmex antillana
Monomorium floricola
Monomorium pharaonis
Odontomachus bauri
Paratrechina longicornis
Pheidole subarmata
Pseudomyrmex subater
Solenopsis geminata
Tapinoma litorale
Tapinoma melanocephalum
Tetramorium bicarinatum
Tetramorium caldarium
Tetramorium lanuginosum
data.frame(`Número de sitios donde fue reportada la especie` = sort(colSums(mc), decreasing = T),
           check.names = F) %>%
  rownames_to_column('Especie') %>% 
  estilo_kable(
    titulo = 'Número de sitios en los que está presente cada especie (orden descendente por número de sitios)', 
    nombres_filas = F, alinear = 'cr')
TABLA 3.3: Número de sitios en los que está presente cada especie (orden descendente por número de sitios)
Especie Número de sitios donde fue reportada la especie
Brachymyrmex heeri 11
Dorymyrmex antillana 11
Pheidole subarmata 9
Paratrechina longicornis 8
Tetramorium lanuginosum 8
Solenopsis geminata 7
Tetramorium bicarinatum 5
Monomorium floricola 4
Cardiocondyla emeryi 3
Tapinoma melanocephalum 3
Brachymyrmex obscurior 2
Cyphomyrmex rimosus 2
Monomorium pharaonis 2
Odontomachus bauri 2
Pseudomyrmex subater 2
Cardiocondyla minutior 1
Cyphomyrmex sp. 1
Tapinoma litorale 1
Tetramorium caldarium 1
data.frame(`Riqueza por sitios` = rowSums(mc),
           check.names = F) %>%  rownames_to_column('Sitio') %>% 
  arrange(desc(`Riqueza por sitios`)) %>% 
  estilo_kable(
    titulo = 'Riqueza por sitios (orden descendente por riqueza)', 
    nombres_filas = F, alinear = 'cr')
TABLA 3.4: Riqueza por sitios (orden descendente por riqueza)
Sitio Riqueza por sitios
H0911 10
H0517 7
H0736 7
H1119 7
H0647 6
H0834 6
H1069 6
H1304 6
H0268 5
H0880 5
H1207 5
H1344 4
H0571 3
H0944 3
H0956 3

La matriz de comunidad analizada se compone de 15 sitios y 19 especies, donde el/los sitio/s más ricos es/son H0911. La/s especie/s más común/es es/son Brachymyrmex heeri y Dorymyrmex antillana y la/s más rara/s es/son Cardiocondyla minutior, Cyphomyrmex sp., Tapinoma litorale y Tetramorium caldarium. El siguiente gráfico de mosaicos muestra la distribución de las especies según sitios.

grafico_mosaico <- crear_grafico_mosaico_de_mc(mc, tam_rotulo = 12) + xlab('Sitios') + ylab('Especie')
grafico_mosaico
Distribución de las especies según sitios

FIGURA 3.1: Distribución de las especies según sitios

3.2 Transformar la matriz de comunidad

Este paso es importante, lo explico aquí

mc_t <- decostand(mc, 'hellinger') #Hellinger, funciona con datos de presencia/ausencia
mc_t %>% estilo_kable(titulo = 'Matriz de comunidad transformada',
                      nombres_filas = T, alinear = 'r')
TABLA 3.5: Matriz de comunidad transformada
Brachymyrmex heeri Brachymyrmex obscurior Cardiocondyla emeryi Cardiocondyla minutior Cyphomyrmex rimosus Cyphomyrmex sp. Dorymyrmex antillana Monomorium floricola Monomorium pharaonis Odontomachus bauri Paratrechina longicornis Pheidole subarmata Pseudomyrmex subater Solenopsis geminata Tapinoma litorale Tapinoma melanocephalum Tetramorium bicarinatum Tetramorium caldarium Tetramorium lanuginosum
H0268 0.45 0.00 0.00 0.00 0.00 0.00 0.45 0.00 0.00 0.00 0.00 0.45 0.45 0.45 0.00 0.00 0.00 0.00 0.00
H0517 0.38 0.00 0.38 0.00 0.00 0.38 0.38 0.38 0.00 0.00 0.38 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.38
H0571 0.58 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.58 0.00 0.00 0.00 0.00 0.58 0.00 0.00 0.00
H0647 0.41 0.00 0.00 0.00 0.00 0.00 0.41 0.41 0.41 0.00 0.00 0.00 0.00 0.41 0.00 0.00 0.41 0.00 0.00
H0736 0.38 0.38 0.00 0.00 0.00 0.00 0.38 0.00 0.38 0.00 0.38 0.38 0.00 0.00 0.00 0.00 0.00 0.00 0.38
H0834 0.41 0.00 0.00 0.00 0.41 0.00 0.00 0.00 0.00 0.41 0.00 0.41 0.00 0.41 0.00 0.00 0.00 0.00 0.41
H0880 0.45 0.00 0.00 0.00 0.00 0.00 0.45 0.00 0.00 0.00 0.00 0.45 0.00 0.45 0.00 0.00 0.00 0.00 0.45
H0911 0.32 0.00 0.32 0.00 0.32 0.00 0.32 0.00 0.00 0.00 0.32 0.32 0.32 0.32 0.00 0.32 0.00 0.00 0.32
H0944 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.58 0.00 0.58 0.00 0.00 0.00 0.00 0.58
H0956 0.00 0.00 0.00 0.00 0.00 0.00 0.58 0.58 0.00 0.00 0.58 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
H1069 0.41 0.00 0.41 0.41 0.00 0.00 0.41 0.00 0.00 0.00 0.00 0.41 0.00 0.00 0.00 0.00 0.41 0.00 0.00
H1119 0.38 0.38 0.00 0.00 0.00 0.00 0.38 0.38 0.00 0.00 0.38 0.00 0.00 0.00 0.00 0.00 0.38 0.00 0.38
H1207 0.45 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.45 0.00 0.45 0.00 0.00 0.45 0.00 0.45 0.00 0.00
H1304 0.00 0.00 0.00 0.00 0.00 0.00 0.41 0.00 0.00 0.00 0.41 0.41 0.00 0.41 0.00 0.00 0.00 0.41 0.41
H1344 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.00 0.00 0.00 0.50 0.00 0.00 0.00 0.00 0.50 0.50 0.00 0.00
# Otras transformaciones posibles con datos de presencia/ausencia
# mc_t <- decostand(mc, 'normalize') #Chord
# mc_t <- decostand(log1p(mc), 'normalize') #Chord
# mc_t <- decostand(mc, 'chi.square') #Chi-square

3.3 Cargar la matriz ambiental

fuente_env_sf <- st_read('data/h3-res-12-no-edificios-3-grupos.gpkg') %>%
  rename(hexagono = indice_propio)
## Reading layer `h3-res-12-no-edificios-3-grupos' from data source 
##   `/home/jose/Documentos/clases_UASD/202302/geo131/repo-gh/manuscrito/data/h3-res-12-no-edificios-3-grupos.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 1490 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 403000 ymin: 2040000 xmax: 404000 ymax: 2040000
## Projected CRS: WGS 84 / UTM zone 19N
env <- fuente_env_sf %>%
  st_drop_geometry() %>% 
  select(-index, -grupo) %>% 
  filter(hexagono %in% rownames(mc)) %>% 
  column_to_rownames('hexagono') %>% 
  mutate(porc_CONS = porc_CONS + porc_EDIF) %>% 
  select(-porc_EDIF)
env %>% estilo_kable(titulo = 'Matriz ambiental', nombres_filas = T, alinear = 'r')
TABLA 3.6: Matriz ambiental
porc_SUEL porc_DOSE porc_CONS
H0268 8.96 47.1 43.92
H0517 43.28 15.9 40.81
H0571 0.00 0.0 100.00
H0647 15.52 42.3 42.15
H0736 25.14 38.6 36.20
H0834 0.00 100.0 0.00
H0880 6.79 37.3 55.91
H0911 0.00 55.8 44.17
H0944 54.75 40.2 5.09
H0956 59.18 0.0 40.82
H1069 7.46 68.7 23.82
H1119 55.13 21.9 22.92
H1207 14.16 58.9 26.99
H1304 0.00 49.0 51.05
H1344 100.00 0.0 0.00

Mapa de hexágonos visitados.

env_sf <- fuente_env_sf %>% 
  select(-index, -grupo) %>% 
  filter(hexagono %in% rownames(mc)) %>% 
  st_transform(4326)
leaflet(env_sf) %>%
  addPolygons(
    popup = ~ paste0(hexagono), stroke = T,
    weight = 2, fillColor = 'red', color = 'black', fillOpacity = 1,
    label = ~ hexagono,
    labelOptions = labelOptions(
      noHide = TRUE, direction = 'auto', offset = c(5, 10), textOnly = T,
      style = list('color' = 'black', 'font-weight' = 'bold',
                   'font-size' = '14px'))) %>%
  addTiles(group = 'OSM') %>%
  addProviderTiles("Esri.NatGeoWorldMap", group="ESRI Mapa") %>%
  addProviderTiles("Esri.WorldImagery", group="ESRI Imagen") %>%
  addProviderTiles("CartoDB.Positron", group= "CartoDB") %>%
  addLayersControl(
    baseGroups = c("CartoDB", "ESRI Imagen", "OSM", "ESRI Mapa"),
    position = 'bottomright',
    options = layersControlOptions(collapsed = FALSE)) %>% 
  addFullscreenControl()

La matriz ambiental se compone de 3 variables de tipo numérico, conteniendo el valor de cada variable para cada uno de los 15 sitios. La siguiente tabla y el gráfico muestran un resumen de los estadísticos básicos de la matriz ambiental.

estad_basicos <- env %>%
  pivot_longer(everything(), names_to = "Variable", values_to = "Valor") %>%
  group_by(Variable) %>%
  summarise(
    Media = mean(Valor),
    Mediana = median(Valor),
    `Desv. Estándar` = sd(Valor),
    Varianza = var(Valor),
    `Error Estándar` = sd(Valor) / sqrt(length(Valor)))
estad_basicos %>% estilo_kable(titulo = 'Matriz ambiental', nombres_filas = F, alinear = 'crrrr')
TABLA 3.7: Matriz ambiental
Variable Media Mediana Desv. Estándar Varianza Error Estándar
porc_CONS 35.6 40.8 25.2 634 6.50
porc_DOSE 38.4 40.2 27.8 775 7.19
porc_SUEL 26.0 14.2 29.9 893 7.72
env %>%
  pivot_longer(everything(), names_to = 'Variable', values_to = 'Valor') %>% 
  group_by(Variable) %>% 
  ggplot() +
  aes(x = Variable, y = Valor, color = Variable, fill = Variable) + 
  # geom_boxplot(lwd = 0.2) + 
  geom_violin(alpha = 0.2, width = 0.8, color = "transparent") +
  geom_jitter(alpha = 0.6, size = 2, height = 0, width = 0.1) +
  geom_boxplot(alpha = 0, width = 0.3, color = "#808080") +
        scale_fill_brewer(palette = 'Set1') +
        theme_bw() +
        theme(legend.position="none")

Las medias calculadas de las variables porc_SUEL, porc_DOSE y porc_CONS son, respectivamente, las siguientes: 26.03, 38.38 y 35.59. La variable que con la media más alta fue porc_DOSE (38.38), y la más baja la obtuvo la variable porc_SUEL (26.03). Por otra parte, la mitad de los sitios midieron menos de 14.16, 40.16 y 40.81, para cada una de las variables porc_CONS, porc_DOSE y porc_SUEL, respectivamente. Finlamente, la variable con mayor dispersión fue porc_SUEL y la de menor dispersión fue porc_CONS.

Una verificación importante que debe realizarse es si las matrices de comunidad y ambiental tienen el mismo numero de filas y si las filas se encuentran en el mismo orden (e.g. consistencia entre matrices, donde cada fila en la matriz de comunidad se refiere al mismo sitio en la ambiental, y viceversa). Esto se puede comprobar por medio de los nombres de columnas y, en este caso, tras realizar la correspondiente comprobación, esta condición se cumple, por lo que podemos continuar adelante con los siguientes análisis

A continuación, realizaré análisis de agrupamiento, ordenación y diversidad, basándome en las indicaciones de Borcard, Gillet, y Legendre (2018), reaprovechando el código contenido en Martínez Batlle (2020).

4 Análisis de agrupamiento

A continuación, el análisis de agrupamiento propiamente. La parte más importante es generar un árbol, a partir de una matriz de distancias, que haga sentido desde el punto de vista de la comunidad y la distribución de las especies. Primero cargaré paquetes específicos de esta técnica y generaré la matriz de distancias.

mc_d <- vegdist(mc_t, "euc")

4.1 Generación de árboles

A continuación, generaré árboles usando distintos métodos. Explico detalladamente estas técnicas en el repo, y en los vídeos (13 a 16) de la lista mencionada arriba “Ecología Numérica con R” de mi canal.

lista_cl <- list(
        cl_single = hclust(mc_d, method = 'single'),
        cl_complete = hclust(mc_d, method = 'complete'),
        cl_upgma = hclust(mc_d, method = 'average'),
        cl_ward = hclust(mc_d, method = 'ward.D2')
)
par(mfrow = c(2,2))
invisible(map(names(lista_cl), function(x) plot(lista_cl[[x]], main = paste0(x, '\n(árbol de evaluación)'), hang = -1)))

par(mfrow = c(1,1))

A continuación, calcularé la distancia y la correlación cofenéticas; esta última, la correlación cofenética,se utiliza como criterio flexible para elegir el método de agrupamiento idóneo, pero no debe usarse de manera estricta. Se supone que el método con la mayor correlación cofenética explica mejor el agrupamiento de la comunidad. Si quieres comprender mejor esta técnica, consulta el vídeo que te referí en el párrafo anterior, así como los libros de referencia. Normalmente, el método UPGMA obtiene la mayor correlación cofenética, pero esto se debe a que su procedimiento de obtención maximiza precisamente dicha métrica. No es recomendable conservar un único método de agrupamiento, normalmente es bueno usar al menos dos. Ward es muchas veces recomendado como método de contraste, por basarse en procedimientos de cálculo muy distintos a los de UPGMA.

map_df(lista_cl, function(x) {
        coph_d <- cophenetic(x)
        corr <- cor(mc_d, coph_d)
        return(corr)
}) %>% t() %>% as.data.frame() %>%
  rownames_to_column %>%
  mutate(rowname = gsub('cl_', '', rowname)) %>% 
  setNames(c('Método de agrupamiento', 'Correlación cofenética')) %>%
  estilo_kable()
TABLA 4.1:
Método de agrupamiento Correlación cofenética
single 0.57
complete 0.69
upgma 0.70
ward 0.66

4.2 Anchura de siluetas

Ahora, calcularé las anchuras de silueta, una métrica que ayuda a determinar en cuántos grupos se organiza la comunidad; las anchuras de silueta no deben usarse como método estricto, y sólo debe usarse de forma flexible para informarnos sobre el número máximo de grupos posibles. Considera las siguientes reglas:

  • El número ideal es 3 grupos, de 4 a 5 grupos es aceptable, 6 o más grupos se considera difícil de interpretar, o es un resultado poco útil; 1 grupo es un resultado sin sentido.
  • Si obtienes distintos grupos, pero uno o varios están compuestos por un único sitio, observa qué ocurre en ese sitio, pues es probable que contenga especie raras sólo presentes en él. En este caso, es recomendable explorar dos alternativas para evitar el grupo formado por un único sitio: ver qué ocurre usando distintos métodos o elegir cortar el árbol en un número de grupos menor.

4.2.1 Anchuras de siluetas para método UPGMA

# UPGMA
anch_sil_upgma <- calcular_anchuras_siluetas(
        mc_orig = mc, 
        distancias = mc_d, 
        cluster = lista_cl$cl_upgma)
u_dend_reord <- reorder.hclust(lista_cl$cl_upgma, mc_d)
plot(u_dend_reord, hang = -1, main = 'Método UPGMA\n(árbol de evaluación)')
rect.hclust(
        tree = u_dend_reord,
        k = anch_sil_upgma$n_grupos_optimo)

resultado_evaluacion_upgma <- evaluar_arbol(u_dend_reord, anch_sil_upgma$n_grupos_optimo)

Tras cortar el árbol, la evaluación practicada concluyó lo siguiente: “Árbol útil para análisis posteriores, siempre que se corte en 2 grupos”

4.2.2 Anchuras de siluetas para método Ward

# Ward
anch_sil_ward <- calcular_anchuras_siluetas(
        mc_orig = mc, 
        distancias = mc_d, 
        cluster = lista_cl$cl_ward)
w_dend_reord <- reorder.hclust(lista_cl$cl_ward, mc_d)
plot(w_dend_reord, hang = -1, main = 'Método Ward\n(árbol de evaluación)')
rect.hclust(
        tree = w_dend_reord,
        k = anch_sil_ward$n_grupos_optimo)

resultado_evaluacion_ward <- evaluar_arbol(w_dend_reord, anch_sil_ward$n_grupos_optimo)

Tras cortar el árbol, la evaluación practicada concluyó lo siguiente: “Árbol útil para análisis posteriores, siempre que se corte en 2 grupos”.

4.3 Remuestreo por bootstrap multiescalar

Una forma alterna de evaluar árboles consiste en usar el remuestreo por bootstrap multiescalar. No me interesa que profundices en ella, sólo presentártela como técnica probabilística para evaluar árboles generados por métodos determinísticos. La técnica es documentada en Borcard, Gillet, y Legendre (2018), de la cual puedes un resumen en este cuaderno y en este vídeo (minuto 51:33). El remuestreo por bootstrap multiescalar valida la robustez de los análisis de agrupamiento tomando múltiples muestras aleatorias de los datos en diferentes tamaños. Este proceso determina qué grupos son consistentemente identificados como clústeres, generando valores de probabilidad aproximadamente insesgados (AU) que son considerados más fiables que las probabilidades de bootstrap tradicionales (BP). Esta técnica ayuda a identificar y confirmar patrones robustos en los datos.

Lo aplicaré primero al árbol generado por el método UPGMA.

# UPGMA
# if(interactive()) dev.new()
cl_pvclust_upgma <-pvclust(t(mc_t),
                           method.hclust = "average",
                           method.dist = "euc",
                           iseed = 99, # Resultado reproducible
                           parallel = TRUE, quiet = T)
# Añadir los valores de p
plot(cl_pvclust_upgma, hang = -1, main = 'Método UPGMA bootstrap\n(árbol de evaluación)')

# Añadir rectángulos a los grupos significativos
lines(cl_pvclust_upgma)
pvrect(cl_pvclust_upgma, alpha = 0.90, border = 4)

Lo aplicaré también al árbol generado por el método Ward.

# Ward
# if(interactive()) dev.new()
cl_pvclust_ward <- pvclust(t(mc_t),
                           method.hclust = "ward.D2",
                           method.dist = "euc",
                           iseed = 99, # Resultado reproducible
                           parallel = TRUE, quiet = TRUE)
# Añadir los valores de p
plot(cl_pvclust_ward, hang = -1, main = 'Método Ward bootstrap\n(árbol de evaluación)')

# Añadir rectángulos a los grupos significativos
lines(cl_pvclust_ward)
pvrect(cl_pvclust_ward, alpha = 0.91, border = 4)

4.4 Conclusión sobre selección de método de agrupamiento y número de grupos

Basado en lo anterior, elegiré un método de agrupamiento y un número de grupos, y lo exportaré a un archivo que posteriormente podré reaprovechar. La lógica empleada para elegir método de agrupamiento y número de grupos, es la siguiente: si el árbol generado por el método UPGMA no es recomendable (por tener grupos formados 2 o menos elementos), pero Ward sí, se usar el árbol generado por el método Ward y el número de grupos idóneo sugerido por la anchura de silueta. Si UPGMA es recomendable pero Ward no lo es, se usar el árbol generado por el método UPGMA, cortado en el número de grupos sugerido por la anchura de siluetas. Si ambos métodos son recomendables y sugieren el mismo número de grupos, se opta por el arbol generado por el método Ward. Si ambos métodos son recomendables pero sugieren un número diferente de grupos, se elige el método que sugiere menos grupos. Finalmente, si ambos métodos, UPGMA y Ward, resultan ser poco idóneos porque generan grupos muy pequeños (dos o menos elementos), se opta, como último recurso, por elegir el árbol generado por el método Ward cortado en 3 grupos.

grupos_seleccionados <- seleccionar_y_cortar_arbol(
  arbol_upgma = lista_cl$cl_upgma, arbol_ward = lista_cl$cl_ward,
  resultado_evaluacion_upgma = resultado_evaluacion_upgma,
  resultado_evaluacion_ward = resultado_evaluacion_ward)
saveRDS(grupos_seleccionados$resultado,
        paste0(fuentes_manuscrito, 'grupos_seleccionados-', params$estudiante,'.RDS'))

Los árboles generados por los métodos UPGMA y Ward no producen grupos compuestos por dos elementos o menos, y ambos podrían cortarse en el mismo número de grupos. Usamos el árbol generado por el método Ward cortado en 2 grupos . El árbol resultante se muestra a continuación:

# Convierte el hclust en dendrograma
dend <- as.dendrogram(grupos_seleccionados$arbol)

# Corta y colorea el dendrograma en k grupos
dend_colored <- color_branches(dend, k=grupos_seleccionados$k)

# Etiqueta los grupos
labels_colors <- labels_colors(dend_colored)
labels(dend_colored) <- paste0(labels(dend_colored), " (",
                               grupos_seleccionados$resultado[grupos_seleccionados$arbol$order],
                               ")")

# Grafica el dendrograma
# par(mar = c(3, 4, 4, 2) + 0.1) # Ajusta los márgenes
plot(
  dend_colored,
  main=paste(
    'Árbol seleccionado\nMétodo',
    grupos_seleccionados$metodo,
    'cortado en',
    grupos_seleccionados$k, 'grupos'),
  xlab = 'Sitios (grupo de pertenencia)')

4.5 Grupos (clústers), variables ambientales

Apliquemos el análisis de agrupamiento a la matriz ambiental. La clave en este punto es que, si la matriz ambiental presenta patrones parecidos a los de la matriz de comunidad, significa que el agrupamiento utilizado hace sentido entre ambos conjuntos de datos (comunidad y hábitat) de forma consistente. Si ambos conjuntos de datos son consistentes, significa que existe algún grado de asociación, aunque sea sólo una mera asociación estadística.

Agrupar los sitios de muestreo de la matriz ambiental según los grupos previamente definidos.

env_grupos <- env %>%
    rownames_to_column('sitios_de_muestreo') %>% 
    mutate(grupos = as.factor(grupos_seleccionados$resultado)) %>%
    pivot_longer(-c(grupos, sitios_de_muestreo), names_to = "variable", values_to = "valor")

Evaluar efectos entre los grupos (“diferencias significativas”). Se utilizan las pruebas estadísticas ANOVA (evalúa homongeneidad de medias) y Kruskal-Wallis (evalúa homogeneidad de medianas). Las tablas están ordenadas en orden ascendente por la columna p_valor_a, que son los p-valores de la prueba ANOVA.

env_grupos_ak <- env_grupos %>%
  group_by(variable) %>%
  summarise(
    p_valor_a = tryCatch(oneway.test(valor ~ grupos)$p.value, error = function(e) NA),
    p_valor_k = tryCatch(kruskal.test(valor ~ grupos)$p.value, error = function(e) NA)
    ) %>%
  arrange(p_valor_a)
env_grupos_ak %>% estilo_kable(alinear = 'crr')
TABLA 4.2:
variable p_valor_a p_valor_k
porc_DOSE 0.00 0.00
porc_SUEL 0.08 0.05
porc_CONS 0.56 0.91

Explora tus resultados.

env_grupos %>% 
        group_by(variable) %>% 
        ggplot() + aes(x = grupos, y = valor, group = grupos, fill = grupos) + 
        geom_boxplot(lwd = 0.2) + 
        scale_fill_brewer(palette = 'Set1') +
        theme_bw() +
        theme(legend.position="none") +
        facet_wrap(~ variable, scales = 'free_y', ncol = 8)

El objetivo de adjuntarle, a la matriz ambiental, el vector de agrupamiento generado a partir de datos de comunidad, consiste en caracterizar ambientalmente los hábitats de los subgrupos diferenciados según su composición. Observa los resultados de las pruebas estadísticas, de los diagramas de caja, y explora tus resultados:

4.6 Especies con preferencia/fidelidad con grupos (clústers)

Análisis de preferencia/fidelidad de especies con grupos (clusters), mediante el coeficiente de correlación biserial puntual (phi).

set.seed(9999)
phi <- multipatt(
  mc,
  grupos_seleccionados$resultado,
  func = "r.g",
  max.order = 1,
  control = how(nperm = 999))
summary(phi)

 Multilevel pattern analysis
 ---------------------------

 Association function: r.g
 Significance level (alpha): 0.05

 Total number of species: 19
 Selected number of species: 2 
 Number of species associated to 1 group: 2 

 List of species associated to each combination: 

 Group A  #sps.  1 
                   stat p.value    
Pheidole subarmata    1   0.001 ***

 Group B  #sps.  1 
                      stat p.value  
Monomorium floricola 0.707    0.02 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Tabla de especies que presentaron asociación con grupos por medio de phi, usando umbral de significancia (umbral_alfa).

tabla_phi_sign <- phi$sign
tabla_phi_sign_alfa <- tabla_phi_sign[phi$sign$p.value < umbral_alfa, ]
data.frame(
  `Nombre de especie` = rownames(tabla_phi_sign_alfa),
  `P-valor` = tabla_phi_sign_alfa$p.value,
  `Grupo de asociación` = gsub('s\\.', '', names(tabla_phi_sign_alfa)[tabla_phi_sign_alfa$index]),
  check.names = F) %>%
  arrange(`Nombre de especie`) %>% 
  estilo_kable(alinear = 'crr')
TABLA 4.3:
Nombre de especie P-valor Grupo de asociación
Monomorium floricola 0.02 B
Pheidole subarmata 0.00 A

5 Técnicas de ordenación

Me basaré en los scripts que comienzan por to_ de este repo, los cuales explico en los vídeos de “Técnicas de ordenación” de la lista de reproducción “Ecología Numérica con R” de mi canal.

5.1 Ordenación no restringida

5.1.1 PCA aplicado a datos de comunidad transformados

pca_mc_t <- rda(mc_t)
summary(pca_mc_t)

Call:
rda(X = mc_t) 

Partitioning of variance:
              Inertia Proportion
Total           0.575          1
Unconstrained   0.575          1

Eigenvalues, and their contribution to the variance 

Importance of components:
                        PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8    PC9   PC10    PC11    PC12
Eigenvalue            0.164 0.0922 0.0703 0.0513 0.0456 0.0384 0.0328 0.0240 0.0166 0.0122 0.00987 0.00904
Proportion Explained  0.285 0.1603 0.1222 0.0891 0.0793 0.0667 0.0570 0.0418 0.0288 0.0212 0.01716 0.01571
Cumulative Proportion 0.285 0.4450 0.5671 0.6562 0.7355 0.8023 0.8593 0.9011 0.9298 0.9511 0.96821 0.98392
                         PC13    PC14
Eigenvalue            0.00655 0.00270
Proportion Explained  0.01138 0.00469
Cumulative Proportion 0.99531 1.00000

Scaling 2 for species and site scores
* Species are scaled proportional to eigenvalues
* Sites are unscaled: weighted dispersion equal on all dimensions
* General scaling constant of scores:  1.68 


Species scores

                               PC1      PC2      PC3     PC4      PC5     PC6
Brachymyrmex heeri       -0.051640 -0.24713 -0.01612  0.1372 -0.28520 -0.0837
Brachymyrmex obscurior    0.061644  0.03692  0.05124  0.1514  0.00994  0.0476
Cardiocondyla emeryi      0.033748 -0.02686  0.02478 -0.0205 -0.20867  0.1995
Cardiocondyla minutior    0.000355 -0.09156  0.05600 -0.0615 -0.05267  0.1458
Cyphomyrmex rimosus      -0.108114  0.00533 -0.09561  0.0404 -0.06563 -0.0409
Cyphomyrmex sp.           0.048579  0.04092  0.02528  0.0767 -0.08887  0.0436
Dorymyrmex antillana      0.213614  0.15996  0.23390 -0.1995 -0.06450  0.0691
Monomorium floricola      0.268835  0.10393  0.25302  0.1187 -0.01928 -0.1339
Monomorium pharaonis      0.027878 -0.01085  0.12079  0.0424 -0.00733 -0.0975
Odontomachus bauri       -0.129281 -0.18782 -0.04396  0.1341  0.09216 -0.0352
Paratrechina longicornis  0.409853  0.21780 -0.21992  0.0475  0.04048  0.0395
Pheidole subarmata       -0.446938 -0.02052 -0.00805 -0.0587  0.03996  0.1511
Pseudomyrmex subater     -0.072879  0.01212 -0.02817 -0.1529 -0.15174 -0.0519
Solenopsis geminata      -0.391947  0.17259  0.05142 -0.1503  0.01168 -0.1909
Tapinoma litorale        -0.036354 -0.16937 -0.00486  0.0580  0.09065  0.0158
Tapinoma melanocephalum   0.188712 -0.05242 -0.34985 -0.1253 -0.01540 -0.0538
Tetramorium bicarinatum   0.159116 -0.31647  0.14266 -0.0657  0.17251  0.0546
Tetramorium caldarium    -0.028701  0.09847 -0.01249 -0.0411  0.08485  0.0467
Tetramorium lanuginosum  -0.247467  0.30142 -0.01609  0.2214  0.01617  0.1012


Site scores (weighted sums of species scores)

           PC1     PC2     PC3     PC4     PC5     PC6
H0268 -0.36616 -0.0740  0.1798 -0.7443 -0.5369 -0.3930
H0517  0.36481  0.3073  0.1898  0.5762 -0.6673  0.3272
H0571  0.43961 -0.2534 -1.0963  0.1415 -0.3105 -0.4407
H0647  0.16294 -0.2735  0.8044 -0.1833  0.0082 -0.9663
H0736  0.03336  0.2139  0.0382  0.5166 -0.0639  0.3113
H0834 -0.64607 -0.1282 -0.2719  0.5291  0.0105 -0.3541
H0880 -0.46278  0.2104  0.1954 -0.0824 -0.2034 -0.0317
H0911 -0.13631  0.2134 -0.5071 -0.3200 -0.6026  0.0903
H0944 -0.72739  0.4258 -0.0753  0.0341  0.5305  0.0449
H0956  0.68639  0.4616  0.3239 -0.0700  0.2456 -0.2193
H1069  0.00247 -0.6365  0.3893 -0.4279 -0.3662  1.0139
H1119  0.42956  0.0634  0.3466  0.6201  0.1386  0.0461
H1207 -0.23073 -1.0750 -0.0308  0.3681  0.5753  0.1001
H1304 -0.19954  0.6846 -0.0868 -0.2858  0.5899  0.3249
H1344  0.64985 -0.1399 -0.3993 -0.6722  0.6523  0.1465
screeplot(
  pca_mc_t,
  bstick = TRUE,
  npcs = length(pca_mc_t$CA$eig)
)

# Biplot
cleanplot.pca(pca_mc_t, scaling = 1, mar.percent = 0.06, cex.char1 = 0.7)

5.1.2 Análisis de correspondencia (CA)

# Realizar el CA
mc_ca <- cca(mc)

Resumen de análisis de correspondencia.

summary(mc_ca)

Call:
cca(X = mc) 

Partitioning of scaled Chi-square:
              Inertia Proportion
Total             2.3          1
Unconstrained     2.3          1

Eigenvalues, and their contribution to the scaled Chi-square 

Importance of components:
                        CA1   CA2   CA3   CA4    CA5    CA6    CA7    CA8    CA9   CA10   CA11   CA12
Eigenvalue            0.431 0.356 0.288 0.267 0.2179 0.1961 0.1464 0.1123 0.0928 0.0750 0.0488 0.0348
Proportion Explained  0.188 0.155 0.125 0.116 0.0947 0.0852 0.0637 0.0488 0.0404 0.0326 0.0212 0.0151
Cumulative Proportion 0.188 0.342 0.467 0.583 0.6781 0.7633 0.8270 0.8758 0.9161 0.9487 0.9699 0.9851
                        CA13    CA14
Eigenvalue            0.0235 0.01085
Proportion Explained  0.0102 0.00472
Cumulative Proportion 0.9953 1.00000

Scaling 2 for species and site scores
* Species are scaled proportional to eigenvalues
* Sites are unscaled: weighted dispersion equal on all dimensions


Species scores

                             CA1     CA2     CA3     CA4     CA5      CA6
Brachymyrmex heeri        0.1309  0.1816  0.0658 -0.0502  0.0117  0.28121
Brachymyrmex obscurior   -0.6495  0.3140 -1.0459  0.5305 -0.5400  0.51512
Cardiocondyla emeryi     -0.3106  0.1477  0.5769 -1.3210  0.4802 -0.03823
Cardiocondyla minutior    0.0727  0.9721  0.7144 -2.8291 -1.3559 -0.59984
Cyphomyrmex rimosus       0.8708 -0.6752  0.6212  0.1168  0.6538  0.74007
Cyphomyrmex sp.          -1.1104  0.2172 -0.2388 -1.0662  2.5027 -0.23414
Dorymyrmex antillana     -0.3974 -0.0203 -0.0795 -0.1037 -0.1603 -0.07862
Monomorium floricola     -1.0233  0.4880 -0.7496 -0.0036  0.5945  0.14314
Monomorium pharaonis     -0.4753  0.2778 -1.1943  0.3419 -0.8951  0.89521
Odontomachus bauri        2.0671  0.9066 -0.1470  0.5129  0.7874  0.00378
Paratrechina longicornis -0.7267  0.0788  0.2739  0.5026  0.1570 -0.30041
Pheidole subarmata        0.6351 -0.3285 -0.0727 -0.1658 -0.2101 -0.11124
Pseudomyrmex subater      0.3146 -0.9059  0.8307 -0.2366 -0.1594  0.90776
Solenopsis geminata       0.4458 -0.8216 -0.1590  0.0156 -0.1494  0.07562
Tapinoma litorale         2.4984  2.4172 -0.2812  0.7246  0.5610 -0.75327
Tapinoma melanocephalum  -0.5903  0.2407  1.8895  1.0099 -0.2157  0.02857
Tetramorium bicarinatum   0.0575  1.1222  0.0301 -0.0556 -0.5386 -0.15898
Tetramorium caldarium     0.0260 -1.4492 -0.5047  0.5021 -0.4236 -2.88409
Tetramorium lanuginosum   0.0844 -0.5520 -0.3304  0.0524  0.2324 -0.09434


Site scores (weighted averages of species scores)

          CA1     CA2     CA3     CA4     CA5     CA6
H0268  0.5234 -1.0654  0.4063 -0.4057 -0.6125  1.0962
H0517 -1.1104  0.2172 -0.2388 -1.0662  2.5027 -0.2341
H0571 -0.9165  0.4696  2.5790  1.8283 -0.0719  0.0159
H0647 -0.4875  0.5752 -1.2070  0.0902 -0.8696  0.9839
H0736 -0.4631 -0.0196 -1.1817  0.5935 -0.9205  0.8065
H0834  1.6358 -0.6040 -0.0129  0.3011  1.0139  0.7608
H0880  0.4167 -0.8664 -0.3997 -0.1888 -0.2530  0.0741
H0911  0.1059 -0.7464  1.2552 -0.0675  0.2937  0.7193
H0944  0.9005 -1.5951 -0.6504 -0.1223 -0.1944 -0.2209
H0956 -1.6593  0.5121 -0.6423  0.4942  0.9042 -0.4010
H1069  0.0727  0.9721  0.7144 -2.8291 -1.3559 -0.5998
H1119 -0.8359  0.6475 -0.9101  0.4674 -0.1595  0.2238
H1207  2.4984  2.4172 -0.2812  0.7246  0.5610 -0.7533
H1304  0.0260 -1.4492 -0.5047  0.5021 -0.4236 -2.8841
H1344 -0.9602  0.9990  1.8344  1.2689 -0.8689 -0.6495

Gráfico de sedimentación o screeplot.

# Screeplot
screeplot(mc_ca, bstick = TRUE, npcs = length(mc_ca$CA$eig))

Representación del biplot.

# Biplot
plot(mc_ca,
     scaling = 1,
     main = "Análisis de correspondencia, escalamiento 1"
)

5.2 Ordenación restringida con modelización

A continuación, el análisis de ordenación propiamente. La parte más importante es el entrenamiento: la función train del paquete caret, contenida en la función my_train, simplifica la selección de variables. Lo más importante: prueba con todas las variables primero, observa las variables que recomienda el modelo final (print_my_train(mod)) y ensaya varias combinaciones de subconjuntos de variables.

mc_t_ren <- mc_t %>%
  rename_all(~ paste('ESPECIE', .x))
env_spp <- env %>% bind_cols(mc_t_ren)
spp <- paste0('`', grep('^ESPECIE', colnames(env_spp), value = T), '`', collapse = ' + ')
my_formula <- as.formula(paste(spp, '~ .'))
set.seed(1); mod <- my_train(
  formula = my_formula, 
  # preproceso = 'scale',
  data = env_spp,
  num_variables = 2)
print_my_train(mod)
$resumen_variables
Subset selection object
3 Variables  (and intercept)
          Forced in Forced out
porc_SUEL     FALSE      FALSE
porc_DOSE     FALSE      FALSE
porc_CONS     FALSE      FALSE
1 subsets of each size up to 2
Selection Algorithm: 'sequential replacement'
         porc_SUEL porc_DOSE porc_CONS
1  ( 1 ) " "       "*"       " "      
2  ( 1 ) "*"       "*"       " "      

$resultados_nvmax
  nvmax  RMSE Rsquared   MAE RMSESD RsquaredSD  MAESD
1     2 0.506    0.436 0.406 0.0373      0.397 0.0625

$mejor_ajuste
  nvmax
1     2
(covar <- grep(
  pattern = '\\(Intercept\\)',
  x = names(coef(mod$finalModel,unlist(mod$bestTune))),
  invert = T, value = T))
[1] "porc_SUEL" "porc_DOSE"
rda_mc_t <- rda(mc_t_ren %>% rename_all(~ gsub('^ESPECIE ', '', .)) ~ .,
                    env %>% select_at(all_of(gsub('\\`', '', covar))), scale = T)

A continuación, el resumen del análisis de redundancia.

summary(rda_mc_t)

Call:
rda(formula = mc_t_ren %>% rename_all(~gsub("^ESPECIE ", "",      .)) ~ porc_SUEL + porc_DOSE, data = env %>% select_at(all_of(gsub("\\`",      "", covar))), scale = T) 

Partitioning of correlations:
              Inertia Proportion
Total           19.00      1.000
Constrained      3.64      0.191
Unconstrained   15.36      0.809

Eigenvalues, and their contribution to the correlations 

Importance of components:
                       RDA1   RDA2   PC1   PC2   PC3    PC4    PC5    PC6    PC7    PC8    PC9  PC10   PC11
Eigenvalue            2.915 0.7204 2.829 2.389 2.035 1.7653 1.5253 1.3345 1.1742 0.8461 0.5821 0.400 0.2678
Proportion Explained  0.153 0.0379 0.149 0.126 0.107 0.0929 0.0803 0.0702 0.0618 0.0445 0.0306 0.021 0.0141
Cumulative Proportion 0.153 0.1913 0.340 0.466 0.573 0.6660 0.7463 0.8165 0.8783 0.9228 0.9535 0.975 0.9886
                        PC12
Eigenvalue            0.2165
Proportion Explained  0.0114
Cumulative Proportion 1.0000

Accumulated constrained eigenvalues
Importance of components:
                       RDA1  RDA2
Eigenvalue            2.915 0.720
Proportion Explained  0.802 0.198
Cumulative Proportion 0.802 1.000

Scaling 2 for species and site scores
* Species are scaled proportional to eigenvalues
* Sites are unscaled: weighted dispersion equal on all dimensions
* General scaling constant of scores:  4.04 


Species scores

                            RDA1    RDA2     PC1      PC2      PC3     PC4
Brachymyrmex heeri       -0.2741  0.5020 -0.4727 -0.04267 -0.09073  0.2647
Brachymyrmex obscurior    0.1172 -0.1334 -0.0153 -0.13311 -0.59564  0.1955
Cardiocondyla emeryi     -0.1490  0.0469 -0.2093 -0.64591  0.49409  0.1776
Cardiocondyla minutior   -0.2779 -0.0346 -0.3934 -0.50246  0.21886 -0.3580
Cyphomyrmex rimosus      -0.5639 -0.0684  0.1452  0.09904  0.32128  0.3594
Cyphomyrmex sp.           0.2081 -0.0126 -0.0195 -0.39668  0.20034  0.5909
Dorymyrmex antillana      0.3429 -0.0808  0.1268 -0.59788 -0.12030 -0.3376
Monomorium floricola      0.4171 -0.0984 -0.0939 -0.37895 -0.29839  0.2730
Monomorium pharaonis     -0.0335  0.0692 -0.0811 -0.18476 -0.66921 -0.0391
Odontomachus bauri       -0.5329 -0.1619 -0.3524  0.50868 -0.00949  0.2864
Paratrechina longicornis  0.6786  0.1127  0.0917 -0.02678  0.09291  0.1669
Pheidole subarmata       -0.6522 -0.0177  0.2759  0.20117  0.03721 -0.2301
Pseudomyrmex subater     -0.1753  0.1882  0.1839 -0.06613  0.33156 -0.1949
Solenopsis geminata      -0.4315  0.0497  0.6200  0.14280 -0.00417 -0.1876
Tapinoma litorale        -0.1870 -0.0306 -0.5237  0.49055 -0.08397  0.0364
Tapinoma melanocephalum   0.4193  0.1694 -0.1371  0.30496  0.47912 -0.0707
Tetramorium bicarinatum   0.0483 -0.3913 -0.6575 -0.00565 -0.13134 -0.3192
Tetramorium caldarium    -0.1092  0.2016  0.4385  0.01711 -0.07764 -0.2574
Tetramorium lanuginosum  -0.2193 -0.1374  0.6397 -0.02288 -0.13035  0.4105


Site scores (weighted sums of species scores)

         RDA1   RDA2     PC1     PC2    PC3     PC4
H0268 -0.7451  1.957  0.5037 -0.0644  0.446 -1.2363
H0517  1.1017  0.344 -0.0821 -1.6705  0.844  2.4882
H0571  1.1171  2.918 -0.7381  1.3147  0.796  0.3741
H0647  0.3637 -0.344 -0.6634 -0.6817 -1.768 -0.6093
H0736  0.2719  0.174  0.2319 -0.3674 -2.088  0.4243
H0834 -2.1777 -0.391  0.2946  0.7992  0.330  1.5562
H0880 -0.6944  0.508  0.8024  0.1382 -0.337 -0.2831
H0911 -0.6179  1.429  0.5832 -0.3747  1.705  0.3755
H0944 -1.0542 -1.130  1.6183  0.9341  0.203 -0.1407
H0956  1.8111 -1.054  0.1907 -0.3504 -0.152 -0.0803
H1069 -0.6882 -0.567 -1.6567 -2.1159  0.922 -1.5074
H1119  1.0942 -1.513 -0.3195 -0.3965 -1.330  0.6974
H1207 -1.3612 -1.220 -2.2053  2.0658 -0.354  0.1534
H1304 -0.0869  0.316  1.8465  0.0720 -0.327 -1.0837
H1344  1.6660 -1.427 -0.4062  0.6975  1.111 -1.1283


Site constraints (linear combinations of constraining variables)

         RDA1    RDA2     PC1     PC2    PC3     PC4
H0268 -0.3682  0.4972  0.5037 -0.0644  0.446 -1.2363
H0517  0.8765 -0.0531 -0.0821 -1.6705  0.844  2.4882
H0571  1.3478  2.4625 -0.7381  1.3147  0.796  0.3741
H0647 -0.1736  0.3534 -0.6634 -0.6817 -1.768 -0.6093
H0736 -0.0124  0.0318  0.2319 -0.3674 -2.088  0.4243
H0834 -2.3452 -0.8334  0.2946  0.7992  0.330  1.5562
H0880 -0.0113  0.9208  0.8024  0.1382 -0.337 -0.2831
H0911 -0.7140  0.6224  0.5832 -0.3747  1.705  0.3755
H0944  0.0112 -1.3802  1.6183  0.9341  0.203 -0.1407
H0956  1.5063 -0.2604  0.1907 -0.3504 -0.152 -0.0803
H1069 -1.1701 -0.1457 -1.6567 -2.1159  0.922 -1.5074
H1119  0.6848 -0.7975 -0.3195 -0.3965 -1.330  0.6974
H1207 -0.7875 -0.1286 -2.2053  2.0658 -0.354  0.1534
H1304 -0.4600  0.8491  1.8465  0.0720 -0.327 -1.0837
H1344  1.6156 -2.1383 -0.4062  0.6975  1.111 -1.1283


Biplot scores for constraining variables

            RDA1    RDA2 PC1 PC2 PC3 PC4
porc_SUEL  0.666 -0.7461   0   0   0   0
porc_DOSE -0.998 -0.0581   0   0   0   0

La varianza ajustada explicada por el modelo.

RsquareAdj(rda_mc_t)$adj.r.squared
[1] 0.0566

Y el factor de inflación de la varianza.

vif.cca(rda_mc_t)
porc_SUEL porc_DOSE 
     1.63      1.63 

Represento el gráfico triplot.

# Triplot
escalado <- 1
plot(rda_mc_t,
     scaling = escalado,
     display = c("sp", "lc", "cn"),
     main = paste("Triplot de RDA especies ~ variables, escalamiento", escalado)
)
rda_mc_t_sc1 <- scores(rda_mc_t,
         choices = 1:2,
         scaling = escalado,
         display = "sp"
  )
# text(mi_fam_t_rda, "species", col="red", cex=0.8, scaling=escalado)
arrows(0, 0,
       rda_mc_t_sc1[, 1] * 0.9,
       rda_mc_t_sc1[, 2] * 0.9,
       length = 0,
       lty = 1,
       col = "red"
)

6 Análisis de diversidad + análisis de agrupamiento abreviado

Me basaré en los scripts que comienzan por di_ de este repo, los cuales explico en los vídeos de “Análisis de diversidad” (vídeos 19 y 20) de la lista de reproducción “Ecología Numérica con R” de mi canal. Dichos vídeos tienen aplicaciones ligeramente diferentes, pues los datos fuente usados en ellos son de abundancia, mientras que los tuyos son de presencia/ausencia.

6.1 Calcular riqueza (e índices)

La principal desventaja de trabajar con registros de presencia, es que la mayoría de los índices de diversidad alpha fueron diseñados originalmente para calcularse a partir de datos de abundancia. Sin embargo, la riqueza de especies, que es el número \(q=0\) de Hill (\(=N_0\) en las columnas que produce la función alpha_div) es un buen proxy sobre la diversidad, y nos ayudará a comparar sitios.

Además de la columna N0 del objeto que generaré en el bloque siguiente, verás que la función alpha_div genera otras columnas; son índices pensados para datos de abundancia, que en este caso no usaremos, pero los muestro para que tengas una visión completa del análisis de diversidad con índices que podría serte de utilidad en el futuro.

Por otra parte, afortunadamente, los métodos de estimación de riqueza de Chao, y los de diversidad beta (al final de esta sección), aprovechan sustancialmente los registros de presencia/ausencia para realizar estimaciones consistentes y fiables.

Una nota adicional. En el análisis de diversidad, es útil (no imprescindible) disponer de un análisis clúster (agrupamiento) básico. Este te servirá para comparar la riqueza observada y la esperada entre hábitats. Por esta razón, combinamos análisis de diversidad con agrupamiento. Sin embargo, si el análisis de agrupamiento generó grupos de dos o menos elementos, dicha comparación no será realizable.

indices <- alpha_div(mc) %>% 
  mutate(sitio = rownames(.)) %>% 
  relocate(sitio, .before = everything())

El objeto mc es la matriz de comunidad de presecia/ausencia. La función alpha_div es un “envoltorio” generado por mí para calcular múltiples índices de diversidad y estimaciones, basada en las funciones de los paquetes SpadeR y iNEXT. Si usásemos datos de abundancia, los índices que calcula la función “alpha_div” serían útiles, pero con registros de presencia/ausencia, como es nuestro caso, sólo la columna N0 (riqueza) nos aportará algún resultado con sentido.

indices %>% 
  kable(booktabs=T) %>%
  kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
  gsub(' NA |NaN ', '', .) #Lista de especies
sitio N0 H Hb2 N1 N1b2 N2 J E10 E20
H0268 H0268 5 1.61 2.32 5 5 5 1 1 1
H0517 H0517 7 1.95 2.81 7 7 7 1 1 1
H0571 H0571 3 1.10 1.58 3 3 3 1 1 1
H0647 H0647 6 1.79 2.58 6 6 6 1 1 1
H0736 H0736 7 1.95 2.81 7 7 7 1 1 1
H0834 H0834 6 1.79 2.58 6 6 6 1 1 1
H0880 H0880 5 1.61 2.32 5 5 5 1 1 1
H0911 H0911 10 2.30 3.32 10 10 10 1 1 1
H0944 H0944 3 1.10 1.58 3 3 3 1 1 1
H0956 H0956 3 1.10 1.58 3 3 3 1 1 1
H1069 H1069 6 1.79 2.58 6 6 6 1 1 1
H1119 H1119 7 1.95 2.81 7 7 7 1 1 1
H1207 H1207 5 1.61 2.32 5 5 5 1 1 1
H1304 H1304 6 1.79 2.58 6 6 6 1 1 1
H1344 H1344 4 1.39 2.00 4 4 4 1 1 1

Los sitios ordenados en función de su riqueza:

indices %>%
  arrange(desc(N0)) %>% 
  kable(booktabs=T) %>%
  kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
  gsub(' NA |NaN ', '', .) #Lista de especies
sitio N0 H Hb2 N1 N1b2 N2 J E10 E20
H0911 H0911 10 2.30 3.32 10 10 10 1 1 1
H0517 H0517 7 1.95 2.81 7 7 7 1 1 1
H0736 H0736 7 1.95 2.81 7 7 7 1 1 1
H1119 H1119 7 1.95 2.81 7 7 7 1 1 1
H0647 H0647 6 1.79 2.58 6 6 6 1 1 1
H0834 H0834 6 1.79 2.58 6 6 6 1 1 1
H1069 H1069 6 1.79 2.58 6 6 6 1 1 1
H1304 H1304 6 1.79 2.58 6 6 6 1 1 1
H0268 H0268 5 1.61 2.32 5 5 5 1 1 1
H0880 H0880 5 1.61 2.32 5 5 5 1 1 1
H1207 H1207 5 1.61 2.32 5 5 5 1 1 1
H1344 H1344 4 1.39 2.00 4 4 4 1 1 1
H0571 H0571 3 1.10 1.58 3 3 3 1 1 1
H0944 H0944 3 1.10 1.58 3 3 3 1 1 1
H0956 H0956 3 1.10 1.58 3 3 3 1 1 1

6.2 Evaluar correlación entre riqueza y variables ambientales mediante matriz de correlación.

En el bloque siguiente, represento gráficamente la correlación entre la riqueza y las variables ambientales mediante un panel de gráficos, que suele llamarse también “matriz de correlación”, expresada gráficamente. Si usases índices de diversidad, como el de Shannon o los números de Hill, también deberías incluirlos en el gráfico; nota que en este ejemplo, sólo uso la riqueza (la función select(N0) se encarga de conservar sólo la riqueza). Esto es lo que debes saber sobre el panel:

  • Presta atención a la primera columna y la primera fila de la matriz, que muestra cómo se correlaciona N0 con las variables ambientales que elijas.

  • La diagonal contiene gráficos de línea que muestra la densidad de la variable en cuestión.

  • Los gráficos del “triángulo superior”, y que contienen el patrón Corr: ####, muestran el valor del coeficiente de correlación de Pearson (\(r\)) entre las variables intersectadas. Si existe un \(|r|\) elevado (es decir, si es muy cercano a -1 o a 1) y la prueba de producto-momento es significativa (si hay uno o varios asteriscos, o un punto, lo es), entonces toma nota de que dicha variable se asocia estadísticamente con la riqueza. Si \(r\) es negativo, la relación es inversa (cuando aumenta la variable, disminuye la riqueza, y viceversa); si es positivo, la relación es directa (cuando aumenta la variable, aumenta también la riqueza).

  • En el “triángulo inferior”, que es un espejo del superior, se sitúan los gráficos de dispersión de las variables intersectadas. Si los puntos siguen un patrón de distribución formando una elipse imaginaria (organizados en torno a una línea recta imaginaria inclinada), entonces existe correlación.

bind_cols(indices %>% select(N0), env %>%
            rename_with(.fn = ~ paste0('AMB_', .))) %>%
  ggpairs(
    labeller = label_wrap_gen(width=10),
    upper = list(continuous = wrap("cor", size = 3))) +
  theme(text = element_text(size = 10))

6.3 “Completitud de muestra” y curva de acumulación

“Completitud”, en porcentajes, según distintos estimadores. Con un 80% de completitud, se considera en general una muestra representativa. Sin embargo, este umbral de 80% no debe tomarse de forma estricta. Sobre todo porque existen métodos refinados que mejoran las estimaciones

riqueza_estimaciones <- data.frame(specpool(mc) %>% select(-matches('.se$'))) %>% 
  select(`Riqueza observada` = Species,
         `Número de sitios` = n,
         `Estimación por Chao (clásico)` = chao,
         `Estimación por jackknife de primer orden` = jack1,
         `Estimación por jackknife de segundo orden` = jack2,
         `Estimación por bootstrap` = boot) %>% 
  pivot_longer(cols = everything(), names_to = 'Variable', values_to = 'Valor') %>%
  mutate(`Cobertura (%)` = Valor / (filter(., Variable == "Riqueza observada") %>% pull(Valor)) * 100) %>% 
  mutate(`Cobertura (%)` = ifelse(Variable %in% c('Riqueza observada', 'Número de sitios'), NA, `Cobertura (%)`))
riqueza_estimaciones %>% estilo_kable(alinear = 'lrr')
TABLA 6.1:
Variable Valor Cobertura (%)
Riqueza observada 19.0
Número de sitios 15.0
Estimación por Chao (clásico) 20.5 108
Estimación por jackknife de primer orden 22.7 120
Estimación por jackknife de segundo orden 22.2 117
Estimación por bootstrap 21.1 111
# POSIBLE ERROR POR BUG NO RESUELTO. PODRÍA APARECER EL SIGUIENTE ERROR:
# Error in if (var_mle > 0) { : valor ausente donde TRUE/FALSE es necesario
# Si te aparece dicho error, entiendo que el problema está relacionado con
# el número de doubletons (quizá la matriz no tiene), dentro de una función
# interna (SpecInciHomo) del paquete SpadeR. Añadir que la versión de SpadeR usada
# en la aplicación Shiny https://chao.shinyapps.io/SpadeR/, no es la misma que 
# la que se encuentra en GitHub ni en el CRAN, pues esa no tiene bug.
df_spader <- data.frame(V1 = as.integer(c(nrow(mc), colSums(mc))))
df_spader %>% estilo_kable(
  titulo = 'Matriz de datos tal como la requiere el paquete SpadeR',
  nombres_filas = T, alinear = 'r', cubre_anchura = T)
TABLA 6.2: Matriz de datos tal como la requiere el paquete SpadeR
V1
1 15
2 11
3 2
4 3
5 1
6 2
7 1
8 11
9 4
10 2
11 2
12 8
13 9
14 2
15 7
16 1
17 3
18 5
19 1
20 8
tryCatch(
  expr = ChaoSpecies(df_spader, datatype = 'incidence_freq',
            k = min(df_spader$V1), conf=0.95),
  error = function(cond) {
      message("Se detectó un error", appendLF = TRUE)
      message('Esta información podría ayudar a depurar: ', cond, appendLF = TRUE)
      message('\nSaliendo...')},
    warning = function(warn) {
      message("Hubo una advertencia", appendLF = TRUE)
      message('Mostrando la advertencia a continuación: ', warn, appendLF = TRUE)},
    finally = {
      message('Estimación de riqueza generada satisfactoriamente')})
## 
## (1) BASIC DATA INFORMATION:
## 
##                                          Variable Value
##     Number of observed species                  D    19
##     Number of sampling units                    T    15
##     Total number of incidences                  U    83
##     Coverage estimate for entire dataset        C 0.959
##     CV for entire dataset                      CV 0.735
## 
##                                                       Variable Value
##     Cut-off point                                            k     1
##     Total number of incidences in infrequent group    U_infreq     4
##     Number of observed species for infrequent group   D_infreq     4
##     Estimated sample coverage for infrequent group    C_infreq 0.152
##     Estimated CV for infrequent group in ICE         CV_infreq     0
##     Estimated CV1 for infrequent group in ICE-1     CV1_infreq     0
##     Number of observed species for frequent group       D_freq    15
## 
##                            Q1
##     Incidence freq. counts  4
## 
## 
## (2) SPECIES RICHNESS ESTIMATORS TABLE:
## 
##                               Estimate   s.e. 95%Lower 95%Upper
##     Homogeneous Model             41.4 155.58     19.5   1091.8
##     Chao2 (Chao, 1987)            20.5   2.04     19.2     30.2
##     Chao2-bc                      19.9   1.50     19.1     27.5
##     iChao2 (Chiu et al. 2014)     20.5   2.04     19.2     30.2
##     ICE (Lee & Chao, 1994)        41.4 155.58     19.5   1091.8
##     ICE-1 (Lee & Chao, 1994)      41.4 155.58     19.5   1091.8
##     1st order jackknife           22.7   2.69     20.1     32.2
##     2nd order jackknife           22.2   4.40     19.4     43.2
## 
## 
## (3) DESCRIPTION OF ESTIMATORS/MODELS:
## 
## Homogeneous Model: This model assumes that all species have the same incidence or detection probabilities. See Eq. (3.2) of Lee and Chao (1994) or Eq. (12a) in Chao and Chiu (2016b).
## 
## Chao2 (Chao, 1987): This approach uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Chao (1987) or Eq. (11a) in Chao and Chiu (2016b).
##      
## Chao2-bc: A bias-corrected form for the Chao2 estimator; see Chao (2005).
##   
## iChao2: An improved Chao2 estimator; see Chiu et al. (2014).
## 
## ICE (Incidence-based Coverage Estimator): A non-parametric estimator originally proposed by Lee and Chao (1994) in the context of capture-recapture data analysis. The observed species are separated as frequent and infrequent species groups; only data in the infrequent group are used to estimate the number of undetected species. The estimated CV for species in the infrequent group characterizes the degree of heterogeneity among species incidence probabilities. See Eq. (12b) of Chao and Chiu (2016b), which is an improved version of Eq. (3.18) in Lee and Chao (1994). This model is also called Model(h) in capture-recapture literature where h denotes "heterogeneity".
## 
## ICE-1: A modified ICE for highly-heterogeneous cases.
## 
## 1st order jackknife: It uses the frequency of uniques to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 2nd order jackknife: It uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 95% Confidence interval: A log-transformation is used for all estimators so that the lower bound of the resulting interval is at least the number of observed species. See Chao (1987).
# Si apareciera ...
# "Error in if (var_mle > 0) { : valor ausente donde TRUE/FALSE es necesario
# ... entonces usar tabla básica de estimación: "riqueza_estimaciones"

Graficaré la curva de acumulación de especies.

mc_general <- mc %>%
  summarise_all(sum) %>%
  mutate(N = nrow(mc)) %>%
  relocate(N, .before = 1) %>%
  data.frame
nasin_raref <- iNEXT::iNEXT(
  x = t(mc_general),
  q=0,
  knots = 2000,
  datatype = 'incidence_freq')
acumulacion_especies <- iNEXT::ggiNEXT(nasin_raref, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0, 80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
acumulacion_especies

Ahora según los grupos previamente seleccionados en el análisis de agrupamiento.

grupos_seleccionados <- readRDS(paste0(
  fuentes_manuscrito, 'grupos_seleccionados-',
  params$estudiante, '.RDS'))
mc_grupos <- mc %>%
  mutate(g = grupos_seleccionados) %>%
  group_by(g) %>%
  summarise_all(sum) %>%
  select(-g) %>% 
  mutate(N = nrow(mc)) %>% 
  relocate(N, .before = 1) %>% 
  data.frame
nasin_raref_general <- iNEXT::iNEXT(
  x = t(mc_grupos),
  q=0,
  knots = 400,
  datatype = 'incidence_freq')
acumulacion_especies_grupos <- iNEXT::ggiNEXT(nasin_raref_general, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0, 80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
acumulacion_especies_grupos

6.4 Contribución de especies a la diversidad beta (SCBD, species contribution to beta diversity) y contribución local a la diversidad beta (LCBD local contribution to beta diversity)

determinar_contrib_local_y_especie(
    mc = mc,
    alpha = 0.05,
    nperm = 9999,
    metodo = 'sorensen')
## $betadiv
## $beta
## SStotal BDtotal 
##   4.115   0.294 
## 
## $SCBD
## [1] NA
## 
## $LCBD
##  H0268  H0517  H0571  H0647  H0736  H0834  H0880  H0911  H0944  H0956  H1069  H1119  H1207  H1304  H1344 
## 0.0559 0.0608 0.0920 0.0658 0.0478 0.0698 0.0378 0.0462 0.0866 0.0885 0.0669 0.0523 0.0912 0.0546 0.0838 
## 
## $p.LCBD
## H0268 H0517 H0571 H0647 H0736 H0834 H0880 H0911 H0944 H0956 H1069 H1119 H1207 H1304 H1344 
## 0.669 0.561 0.110 0.446 0.812 0.383 0.953 0.856 0.166 0.143 0.452 0.757 0.117 0.704 0.191 
## 
## $p.adj
## H0268 H0517 H0571 H0647 H0736 H0834 H0880 H0911 H0944 H0956 H1069 H1119 H1207 H1304 H1344 
##     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1 
## 
## $method
## [1] "sorensen"     "sqrt.D=FALSE"
## 
## $note
## [1] "Info -- D is Euclidean because beta.div outputs D[jk] = sqrt(1-S[jk])"
## [2] "For this D functions, use beta.div with option sqrt.D=FALSE"          
## 
## $D
## [1] NA
## 
## attr(,"class")
## [1] "beta.div"
## 
## $especies_contribuyen_betadiv
## [1] NA
## 
## $sitios_contribuyen_betadiv
## character(0)
## 
## $valor_de_ajustado_lcbd
## H0268 H0517 H0571 H0647 H0736 H0834 H0880 H0911 H0944 H0956 H1069 H1119 H1207 H1304 H1344 
##     1     1     1     1     1     1     1     1     1     1     1     1     1     1     1 
## 
## $sitios_contribuyen_betadiv_ajustado
## character(0)

Referencias

Borcard, Daniel, François Gillet, y Pierre Legendre. 2018. Numerical ecology with R. Springer.
Martínez Batlle, José Ramón. 2020. biogeografia-master/scripts-de-analisis-BCI: Long coding sessions (versión v0.0.0.9000). Zenodo. https://doi.org/10.5281/zenodo.4402362.